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Abstract 

Stable isotopes ((5 15 N and <5 18 0) of the greenhouse gas N 2 0 provide information about the sources and processes leading to 
N 2 0 production and emission from aquatic ecosystems to the atmosphere. In turn, this describes the fate of nitrogen in the 
aquatic environment since N 2 0 is an obligate intermediate of denitrification and can be a by-product of nitrification. 
However, due to exchange with the atmosphere, the 3 values at typical concentrations in aquatic ecosystems differ 
significantly from both the source of N 2 0 and the N 2 0 emitted to the atmosphere. A dynamic model, SIDNO, was developed 
to explore the relationship between the isotopic ratios of N 2 0, N 2 0 source, and the emitted N 2 0. If the N 2 0 production rate 
or isotopic ratios vary, then the N 2 0 concentration and isotopic ratios may vary or be constant, not necessarily 
concomitantly, depending on the synchronicity of production rate and source isotopic ratios. Thus prima facie 
interpretation of patterns in dissolved N 2 0 concentrations and isotopic ratios is difficult. The dynamic model may be used to 
correctly interpret diel field data and allows for the estimation of the gas exchange coefficient, N 2 0 production rate, and the 
production-weighted 3 values of the N 2 0 source in aquatic ecosystems. Combining field data with these modelling efforts 
allows this critical piece of nitrogen cycling and N 2 0 flux to the atmosphere to be assessed. 
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Introduction 

Nitrous oxide (N 2 0) is a powerful greenhouse gas, 298 times 
more potent than CO z over a 100-year time line [1]. Atmospheric 
N 2 0 concentrations have been increasing at a rate of 0.25%/year 
over the last 150 years [2]. Consequently, the global N 2 0 budget 
has been the subject of intensive research efforts over the past few 
decades. N 2 0 is produced through multiple microbial pathways: 
hydroxylamine oxidation during nitrification and as an obligate 
intermediate during denitrification and nitrifier-denitrification. 
Because these pathways of N 2 0 production have different stable 
isotopic enrichment factors, isotopic analysis of N 2 0 can 
potentially distinguish N 2 0 produced through different pathways 
or from different sources [3] . Identifying N 2 0 sources will provide 
insights on the fate of N at the ecosystem-scale (e.g., [4—6]). The 
isotopic ratios of N 2 0 produced in soil environments (e.g., [7-1 1]), 
and in aquatic environments (e.g., [12-18]) have been measured to 
some extent. Although N 2 0 production in rivers and estuaries is a 
significant portion of the global N 2 0 budget (approximately 
1.5 TgN/year, [19]), few studies report isotopic data for rivers 
[5,20,21]. 

In ice-free aquatic ecosystems, the <5 15 N and <5 18 0 of dissolved 
N 2 0 is affected by gas exchange with the atmosphere. As a result, 



the isotopic ratios of dissolved N 2 0 are not equal to those of the 
N 2 0 produced within the aquatic ecosystem and continue to 
change as atmospheric exchange (both ingassing and outgassing) 
occurs. In addition, isotopic fractionation during influx and efflux 
causes the isotopic ratios of N 2 0 flux emitted to the atmosphere to 
be different than that of the dissolved N 2 0 [22]. Thus, the simple 
method of calculating the instantaneous isotopic ratios of the N z O 
flux by taking measured dissolved isotopic ratios, adding an 
equilibrium isotope fractionation, and applying them to measured 
flux rates is inappropriate. Adjustments of measured isotopic ratios 
are necessary to understand the isotopic ratios of both produced 
and emitted N 2 0. 

In this paper, we present a dynamic model of the stable isotopic 
composition of both the dissolved and emitted N 2 0 in aqueous 
systems. We apply this model to two different measured diel 
patterns of the isotopic ratios of N 2 0 in an aquatic ecosystem. We 
use the model to elucidate the relationship between the isotopic 
ratios of source, dissolved, and emitted N 2 0, to allow for improved 
interpretation of dissolved N 2 0 isotope data. Ultimately, a 
process-based understanding on N cycling with aquatic ecosystems 
may be developed based on interpretation of N cycling processes. 
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Materials and Methods 

Stable Isotopes of N 2 0 

N z O is an asymmetric molecule: the most abundant isotopo- 
logues of N 2 0 are [ 14 N 14 N 16 0], [ 15 N 14 N 16 N], [ 14 N 15 N 16 0] and 
[ 14 N 14 N ls O]. The isotopic ratios, 15 N: 14 N and 18 0: lfi O, are: 



5 R = 



[ l5 N 14 N 16 0] + [ 14 N I5 N 16 0] _ [ 15 N 2 0] 



2[ 14 N I4 N 16 0] 



4 N 2 0] 



(1) 



18r J' 4 N 14 N 18 0]_[N 2 18 0] 



rl4 N 14 N 16 0] [N 2 lb O] 



(2) 



where [ 14 N 14 N 16 0], [ 15 N 14 N 16 0], [ 14 N 15 N 16 0] and [ 14 N 14 N ls O] 
represent the concentrations of the various N 2 0 isotopologues. 
Note that 15R is the bulk 15 N: 14 N ratio and represents an average 
ratio of the two 15 N isotopomers and isotopic ratios are reported as 
8 'N relative to air and <5 ls O relative to VSMOW. Although the 
isotopic ratio of the 15 N isotopomers can be measured (e.g., [23— 
25]), the gas exchange fractionation factors are not affected by the 
intramolecular distribution of 5 N [22] . Many laboratories cannot 
measure the intramolecular distribution of 15 N and analysis of the 
bulk L, N: 14 N ratio of N z O is more common [26]. Here, we 
confine our analysis to bulk 15 N: 14 N ratios and use 15 N 2 0 to 
represent the average abundance of the two 15 N isotopomers. The 
same approach could easily be extended to consider each 
isotopologue separately. 

Dynamic Isotope Model for Dissolved N 2 0 

A simple three box model (SIDNO, Stable Isotopes of Dissolved 
Nitrous Oxide) was created using Stella modelling software 
(version 9.1.4, http://www.iseesystems.com) in order to study 
the relationships between the isotopic ratios of source, dissolved 
and emitted N 2 0 (model file is available at https://github.com/ 
jjvenky /SIDNO and by contacting the corresponding author). 
This model is an adaptation of the isotopic gas exchange portion of 
the PoRGy model [27], which successfully modelled diel isotopic 
ratios of 0 2 resulting from photosynthesis, respiration, and gas 
exchange in aquatic ecosystems. One key difference is photosyn- 
thetically produced 0 2 in PoRGy has a S O value fixed by the 
H 2 0 molecules, whereas SIDNO has N 2 0 production <5 15 N and 
<5 1!i O values that can vary independently of each other and of N 2 0 
production rate in order to simulate variability in nitrification and 
denitrification. 

One box in SIDNO is used for the total mass of dissolved N 2 0 
and two additional boxes for the dissolved masses of the two heavy 
isotopologues ( L, N 2 0 and N 2 18 0). The boxes are open to the 
atmosphere for gas exchange, are depth agnostic, and each box 
can gain N 2 0 via a production term; there is no N 2 0 
consumption term since the 3 values of N 2 0 are largely controlled 
by the production pathways [28,29] though certain waters can 
exhibit significant N 2 0 reduction to N 2 [30,31]. The masses and 
magnitude of the flows of 15 N 2 0 and N 2 18 0 relative to bulk N 2 0 
are used to calculate the isotopic composition of source, dissolved, 
and emitted N 2 0. Although isotopic ratios are used in the model, 
we discuss 5 values that are common for reporting isotopic ratios. 
N 2 0 production rate and its S values are user-defined and can be 
adjusted for diel patterns in N 2 0 production that may be caused 
by variable 0 2 levels [32-37]. 



Stable Isotope Dynamics of Gas Exchange 

The 5 values of the net gas exchange flux are controlled by the 
kinetic fractionation factors for evasion (oc ev = R eva ded I Rdbsohed > 
0.9993 for <5 13 N and 0.9981 for <5 18 0) and invasion 
{a in = R inmded /Rg as , 1.0000 for <S 15 N and 0.9992 for cS 18 0) [22]. 
These two a values are related to the equilibrium fractionation 
factor: a eq = R gas / 'Rdissohed = ««v / '<*-m (0.99925 for <S 15 N and 
0.99894 for <5 18 0) and are independent of temperature over the 
range of 0 C C to 44.5 J C [22]. 

The 5 values of tropospheric N 2 0 are 6.72 %o+ 0.12%o for 
<5 15 N and 44.62%o ± 0.21%o for <5 ls O [38]. Therefore, at 
equilibrium, dissolved N 2 0 has dissolved S values slightly greater 
than these at 7.48%o and 45.73%o, respectively. 

In the model, net N 2 0 flux between the atmosphere and 
dissolved phase was calculated using the thin boundary layer 
approach as: 



Flux = k (pn 2 o x k H - [N 2 0] dissohed j 



(3) 



where the N 2 0 flux is calculated in mol/m 2 /h, k is the user- 
modifiable gas exchange coefficient (m/h), is the partial pressure of 
tropospheric N 2 0 (assumed to be 320 ppbv from data provided by 
the ALE GAGE AGAGE investigators, [39,40]), k H is the Henry 
constant for N 2 0 (mol/atm-m' 5 ), and [N 2 0]^ ssofr( , rf is the dissolved 
concentration of N 2 0 (mol/m' 5 ). ku is a function of water 
temperature [41]: 



= 0.025e 



\T 298.15/ 



(4) 



where T is temperature in kelvins. 

Gas exchange is a two-way process. The net N 2 0 flux rate (the 
difference between the invasion and evasion rates) depends on the 
dissolved N 2 0 concentration. When a solution is at equilibrium 
with the atmosphere, the invasion and evasion rates will be equal, 
and the net flux will be zero. 

As with the bulk N 2 0 flux, the flux of the heavy isotopologues 
( 15 N 2 0 and N 2 18 0) can be calculated by including the kinetic 
fractionation factors for N 2 0 (adapted from [27]): 

Flux 15 K 2 0 = k^ xp 15 ^ Q x k H -a\l[ l5 K 2 0] dbsohed ) (5) 



Flux^O = k $ xp^n Q xk H - all [N 2 18 0] fl 



(6) 



where />i5 N2 o all( f 18 o are t ' le P 31 ^ pressures of N 2 0 and 
N 2 1!i O. 



Results 

Test of Model Performance 

To test the ability of SIDNO to reproduce observed isotopic 
data, input parameters (N 2 0 production rate, N 2 0 S values, and 
k) were set to replicate a series of experiments designed to derive 
fractionation factors for N 2 0 gas exchange [22]. In these 
experiments, degassed water was exposed to N 2 0 gas of known 
isotopic ratios in a sealed container to varying degrees of 
saturation. 
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Figure 1 . Comparing the model output to the experimental data of [22]. The coefficient of determination for experimental data and SIDNO 
model outputs were comparable to those of [22], J? 2 = 0.77 for <5 15 N and R 2 = 0.82 for <5 18 0. Precision of measurements for the experimental data was 
+ 0.05%o for 5 15 N and ± 0.1 %o for S™0. 
doi:1 0.1 371 /journal.pone.0090641 .g001 



Modelled dissolved N 2 0 concentration and <5 values increased 
in response to gas exchange (Figure 1). The model fit to the 
experimental data is comparable to the original best-fit derivations 
(R 2 = 0.77 for <5 13 N and R 2 = 0.82 for 5 18 0 for both the original fit 
[22] and the SIDNO fit) (Figure 1). The initial isotopic 
composition of dissolved N z O was identical to the gas phase 
() 1S N value, but the <5 18 0 of dissolved N 2 0 was slightly less than 
the gas phase <5 O value. Ultimately, at 100% saturation the S 
values of the dissolved N 2 0 were greater than those of the gas 
phase as a result of a eq . The model successfully simulated the 
kinetic and equilibrium fractionations during gas exchange under 
the experimental conditions. 

Next, SIDNO was used to provide insight into the effect of 
degassing on the 3 values of dissolved and emitted N 2 0. Here the 
results of two model runs with the same initial N 2 0 concentration 
but different initial S values of dissolved N 2 0 were compared 
(Figure 2). As N 2 0 saturation declined both the dissolved <5 15 N 
values and instantaneous <5 15 N values of the emitted N 2 0 
remained relatively constant, dissolved S O values and instanta- 



neous S O values of the emitted N 2 0 varied by about 10, when 
the solution was very supersaturated (>300% saturation). The <5 
values rose quickly as the system approached 100% saturation. 
Because the light isotopologue diffuses out of solution faster than 
the heavy isotopologue, the instantaneous S values of the emitted 
N 2 0 were always less than the concomitant 5 values of dissolved 
N 2 0. The isotopologues of N 2 0 reached equilibrium indepen- 
dently of each other and therefore the total mass emitted for each 
isotopologue and rate of change depended on the initial 
concentration and 3 values. The retention of N 2 0 in the dissolved 
phase caused the 5 values of the mass emitted to differ from those 
of total mass production. However, when initial dissolved N z O 
concentrations were high (> 1000% saturation) the S values of the 
total N 2 0 emitted were similar to the S values of dissolved N 2 0 
because the mass of N 2 0 lost is very much larger than the N 2 0 
that remained dissolved. The value of k did not affect the gas 
exchange trajectories only the speed at which the system reached 
equilibrium. 
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Figure 2. <> 15 N and <i 1s O trajectories for dissolved and emitted N 2 0 in two supersaturated solutions with zero N 2 0 production. Initial 
dissolved isotopic values for the two dissolved N 2 0 solutions were (5 15 N = -50% 0 , <5 l8 O = 10% 0 , and c5 15 N = -10% 0 , .5 18 0 = 30%o. Both runs used an 
initial dissolved N 2 0 concentration of 1500% saturation. Note that in the <5 18 0 versus <5 15 N plot, the dissolved N 2 0 curves do not pass through the 
tropospheric N 2 0 value due to the small equilibrium isotope effect. 
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N z O isotope data are often plotted as <5 ls O versus <5 15 N to 
elucidate relationships between the various sources and tropo- 
spheric N 2 0 [38]. The trajectories on these plots (Figure 2C) were 
dictated by the 3 values of the source relative to the constant 
atmospheric value and the a. values. Note that some plots in the 
literature differ due to different reference materials for the <5 ls O 
scale (VSMOW and atmospheric 0 2 ). 

Modelling Scenarios with Steady State Production of N 2 0 

The SIDNO model can be used to probe the stable isotope 
dynamics of N 2 0 in a variety of situations that may be 
encountered in aquatic environments to elucidate the relationship 
between the N 2 0 source (a function of N cycling processes), 
dissolved (the easily measured component), and emitted (of 
consequence for greenhouse gas production and global N and 
N 2 0 cycle). 

In the steady-state production of N 2 0 (constant rate and 3 
values), by definition, the 5 values of N 2 0 production must be the 
same as those of the emitted N 2 0. As a result, the 3 values of the 
dissolved N z O cannot equal that of the source (or emitted) N 2 0 at 
steady state because the dissolved N 2 0 must be offset from the 
emitted N z O by at least the a ev values. As the steady-state 
production rate was increased, the steady-state N 2 0 concentration 
increased and the dissolved 3 values approached but did not equal 
the source (Figure 3). Even at moderate supersaturations 
(< 1000%) the effect of atmospheric N 2 0 equilibration on the 3 
values of dissolved N 2 0 cannot be ignored. 

At steady state, the 3 values of the emitted N 2 0 must be equal 
to the source; the large difference between source/emitted and 
dissolved N 2 0 underscores the importance of adjusting the 
measured 8 values of dissolved N 2 0 in order to determine aquatic 
contributions of N 2 0 to the atmosphere or N 2 0 sources. This is 



critical when using dissolved measurements of N 2 0 to constrain 
the global isotopic N 2 0 budget, but not been done in most studies, 
e.g., [16,42-44] but see [45]. 

Modelling Scenarios with Variable Production of N 2 0 

The relationship between the S values of source, dissolved, and 
emitted N 2 0 are much more complicated when N 2 0 production 
is variable rather than when it is constant. N 2 0 production may 
vary with respect to production rate and/ or 3 values; in many 
aquatic environments, N 2 0 production is not likely to be constant. 
The N 2 0 production processes, nitrification and denitrification, 
are sensitive to redox conditions, which can be highly variable, due 
to diel changes in dissolved 0 2 concentration, flow regime, etc. For 
example, [34] observed diel changes in the denitrification rate in 
the Iroquois River and Sugar Creek (Midwestern USA) and found 
that the denitrification was consistently greater during the day 
than night. The relative importance of nitrification and denitri- 
fication can change in response to the diel oxygen cycle: e.g., [46] 
observed a change from daytime nitrification to nighttime 
denitrification in a subtropical eutrophic stream. Coupling of 
N 2 0 and 0 2 diel cycles has been observed in agricultural and 
waste-water treatment plant (WWTP) impacted rivers [36]. Since 
fractionation factors and substrates are different for nitrification 
and denitrification, ecosystem-scale fractionation factors may be 
rate and process dependent, and the 3 values of N 2 0 production 
in a given ecosystem may not be constant over a diel cycle. 

To simulate the diel variability, various scenarios were modelled 
by adjusting either production rate and/ or the associated 3 values. 
The variabilities in these input parameters were driven by a sine 
function with a 24 h period similar to a dissolve 0 2 curve. In all 
scenarios, the chosen range of production rates was based on 
published N 2 0 flux rates (Table 1) and varied from 1 to 
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scenarios where the S values of source N 2 0 was variable, the sine 
function for the S values was synchronized so that maxima and 
minima <5 L, N values coincided with those of <5 1!i O. This was done 
for simplicity, and because, in general, nitrification yields N 2 0 
with lower <5 L 'N and (5 1h O values than denitrification (e.g., 
[10,48]). Nevertheless, scenarios with greater amounts of N 2 0 
reduction to N 2 can be modelled by increasing the source <5 15 N 
and <5 1!i O values to those appropriate for any given ecosystem. 
Model scenarios were run until the output parameters (i.e., N 2 0 
saturation and the 5 values of dissolved, source, and emitted N 2 0) 
reached dynamic steady state: model output was not constant over 
24 h but the diel patterns on successive days were repeated. 

Model Scenario #1: Variable Production Rate, Constant 
Isotopic Composition of Source 

In scenario #1 (Table 2), the S values of source N 2 0 were held 
constant and the production rate was variable. An example of such 
a system may be N 2 0 production via denitrification in river 
sediments with abundant NO^ . Denitrification rates in rivers have 
been observed to fluctuate in response to the diel 0 2 cycle [34] . If 
the fractionation factors for denitrification are not rate dependent, 
the resulting N 2 0 production rate would be variable but the 
source c5 values of N 2 0 values could be constant. 

Here, the maximum concentration lagged approximately 2.75 h 
behind the maximum N 2 0 production rate, a function of the 
magnitude of the gas exchange coefficient, cf. [49] . The 8 values 
for the instantaneously emitted N 2 0 were relatively constant and 
very similar to the N 2 0 source (within 0.4%o for <5 Ll N and 1.196a 
for (5 1 O, Figure 4, Table 3). However, the 5 values of dissolved 
N 2 0 were more variable, spanning 16%o for (5 15 N and 10%o for 
c> 1!! 0. Thus, a change in the <5 values of dissolved N 2 0 can be 
driven simply by a change in production rate and not necessarily a 
change in the 3 values of the source. Since the system was at 
dynamic steady state, the average 5 values of the emitted N 2 0 
were identical to the average 8 values of the source. This must be 
true in all steady-state cases to conserve the mass of the N 2 0 
isotopologues. 

In some aquatic systems, the N 2 0 production rate may remain 
constant with time but the 5 values of the source may change with 
time. In rivers or lakes without a strong diel 0 2 cycle, sediment 
denitrification may produce N 2 0 at an approximately constant 
rate. Denitrification rate may also be independent of water column 
NOj~ concentration if limited by factors other than diffusion in the 
sediments. The S values of the source N 2 0 may thus change if the 



Table 1. Summary of relevant published data on N 2 0 production in aquatic environments. 





Location 


Range of N 2 0 
Saturation (%) 


Range of N 2 0 
Flux (/imol/m^h 1 ) 


Range of k 
values (m/h) 


Reference 


Ohio River, OH, US 


95 to 745 


5 to 90 




[12] 


5 agricultural streams, ON, CA (over 2-3 years) 


14 to 1700 


-1 to 91.7 


0.002 to 0.59 


[5] 


10 agricultural streams, ON, CA (over 17 diel cycles) 


30 to 2570 


-0.33 to 52.1 


0.004 to 0.30 


[45] 


Bang Nara River, TH 


1 70 to 2000 






[20] 


Lll River, NZ 


201 to 404 


1.35 to 17.9 


0.13 to 0.82 


[58] 


Lll River, NZ 


402 to 644 


0.46 to 0.89 


14.76 


[33] 


Seine River, FR 




2.2 to 5.2 


0.04 to 0.06 


[59] 


Canal Two, Yaqui Valley, MX 


100 to 6000 


0 to 34.9 


0.3 to 0.6 


[46] 


agricultural stream, UK 


100 to 630 


0 to 37.5 




[60] 


Grand River, ON, CA 


38 to 8573 


-1.4 to 173.6 


0.06 to 0.35 


[36,37] 



doi:1 0.1 371 /journal.pone.0090641 .t001 
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5 mol/mVh 1 (Table 2), which was between the diel variation in 
N 2 0 flux observed by [33] and [46]. Temperature was held 
constant at 20°C. The value of k was varied between 0.1 and 
0.3 m/d (Table 2), within the range observed in other river studies 
(Table 1). The combination of production rates and k values were 
chosen to produce N 2 0 between 150% and 500% saturation 
(Table 2) coinciding with the range of published data (Table 1). 
The range of 5 values used for the N 2 0 source (Table 2) was 
within published values from various field studies [47]. For 
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Table 2. Summary of input parameters for the SIDNO model scenarios for non-steady state production of N 2 0. 





Scenario # Results Figure 


k 


N 2 0 Source Production Rate 


N 2 0 Source <5 15 N and <5 18 0 Values 




(m/k) 


(/imol/m^h 1 ) 


Wo) 


Variable Production Rate, Constant isotopic Composition 


of Source 






1 4 


0.3 


1 to 5 


-50, 10 


Constant Production Rate, Variable Isotopic Composition 


of Source 






2 5 


0.3 


3 


-50, 10 to -30, 10 


3 6 


0.1 


3 


-50, 10 to -30, 10 


Variable Production Rate, Variable Isotopic Composition 


of Source 






4* 7 


0.3 


1 to 5 


-50, 10 to -30, 10 


5»* 8 


0.3 


1 to 5 


-50, 10 to -30, 10 


6* 9 


0.1 


1 to 5 


-50, 10 to -30, 10 



*Maximum production rate coincides with the lowest source S values. 
**Maximum production rate coincides with the highest source S values. 
doi:1 0.1 371 /journal.pone.0090641 .t002 
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Figure 4. Model scenario #1 - isotopic composition of dissolved and emitted N 2 0 with a variable production rate and constant 
isotopic composition of the source. Note, in panel D, the data points for emitted N 2 0 are masked by the data point for source N 2 0. 
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Table 3. Summary of SIDNO output for model scenarios simulating non-steady state production of N 2 0. 



Dissolved N 2 0 






Emitted N 2 0 








Scenario # Saturation <5 15 N, S ,e O 


A<5 15 N 


Ad ,a O 


<5 15 N, <5 18 0 




Atj 15 N 


A<5 18 0 


(%) (%o) 


(%>) 


(%.) 


(%>) 




(%o) 


(%o) 


Variable Production Rate, Constant isotopic Composition of Source 


1 153 to 263 -27.8, 24.6 to -12.1, 34.4 


37.9 


24.4 


-49.6, 9.5 to - 


-50.2, 11.1 


0.2 


0.5 


Constant Production Rate, Variable Isotopic Composition of Source 


2 208 -19.6, 29.4 to -3.7, to 37.4 


30.4 


19.4 


-45.3, 12.3 to 


-14.7, 27.7 


4.7 


2.3 


3 423 -26.1, 24.8 to -15.1, 30.3 


23.9 


14.8 


-37.2, 16.4 to 


-22.8, 23.6 


12.8 


6.4 


Variable Production Rate, Variable Isotopic Composition of Source 


4* 1 53 to 263 -25.8, 25.6 to - 1 .2, 39.8 


24.2 


15.6 


-46.9, 11.3 to 


- 1 8.0, 26.5 


8 


3.5 


5** 1 53 to 263 - 1 1 .2, 33.8 to -5.0, 36.1 


38.8 


23.8 


-41.6, 14.6 to 


-13.4, 28.1 


8.4 


4.6 


6* 345 to 501 -31.6, 21.6 to -18.4, 29.3 


18.4 


11.6 


-42.1, 13.7 to 


-29.4, 20.6 


19.4 


9.4 



Temporally variable parameters are given as a range. A<5 15 N and At> 18 0 (A = (5 — S) are the maximum difference between the range of source N 2 0 and the range for the 
model output parameter. 

*Maximum production rate coincides with the lowest source S values. 
**Maximum production rate coincides with the highest source S values. 
doi:1 0.1 371 /journal.pone.0090641 .t003 
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Figure 5. Model scenario #2 - Isotopic composition of dissolved and emitted with a constant production rate and variable isotopic 
composition of the source. 
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Figure 6. Model scenario #3 - Isotopic composition of dissolved and emitted N 2 0 with a constant production rate and variable 
isotopic composition of the source, k is reduced from 0.3 m/h to 0.1 m/h. 
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8 values of the NOj~ substrate changed with time. For example, 
many studies have shown that the 8 values of residual NO^ 
increase during denitrification [50]. Similarly, NO^ from 
WWTPs may have different S values than agricultural runoff 
and diel changes in WWTP release may result in changing 8 
values of NO^~ . Changes in N cycling may also vary on a diel basis 
but result in fortuitously similar N 2 0 production rates due to, for 
example, changes in the N 2 0:N 2 ratio of denitrification or changes 
in the relative importance of nitrification and denitrification. Thus 
changes in 8 values of the N z O source do not necessarily indicate 
changes in N 2 0 production rates. 

Model Scenario #2: Constant Production Rate, Variable 
Isotopic Composition of Source 

In scenario #2, when the N 2 0 production rate was held 
constant and the 8 values of the source varied with time (from 
-50%o to - 10%o for <5 15 N and from 10%o to 30%o for <5 18 0), the 8 
values of the dissolved N 2 0 was also much farther from that of the 
source than the dissolved N 2 0 due to the effects of atmospheric 
exchange and the emitted N 2 0 varies linearly between the two 
source values. In contrast, the dissolved N 2 0 is parallel but offset 
from the line connecting the two sources (Figure 5, Table 3). The 
maximum difference between emitted and source N 2 0 was 4.7%o 



for 8 N and 2.3%o for 8 O. The dissolved and emitted 8 values 
also lagged 2.75 h behind the source as a result of gas exchange (as 
above). Since the system was at dynamic steady state, the average 
8 values of the emitted N 2 0 were identical to the average 8 values 
of the source. 



Model Scenario #3: Constant Production Rate, Variable 
Isotopic Composition of Source 

To examine the effects of varying k on the scenario of constant 
N 2 0 production with variable isotopic signature of the source, k 
was reduced from 0.3 m/h (scenario #2) to 0. 1 m/h (scenario #3; 
Figure 6, Table 3). The 8 values for the emitted N 2 0 were centred 
between the sources N z O values, but dissolved N 2 0 8 values were 
farther from tropospheric N 2 0 than the high-A: scenario #2 
(Figure 6 D). 

The effect of reducing k was an increase in N 2 0 concentration 
with the same production rate and a shift in the 8 values of 
dissolved N 2 0 toward the source values. Reducing k also 
dampened the response between the instantaneous 8 values of 
the emitted N 2 0 and the 8 values of the source. As above, the lag 
time between the 8 values of the source and emitted N 2 0 
increased as k decreased. The total range of the source and 
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emitted S values decreased. The difference between the source and 
emitted <5 values was 12.8%o for <5 I5 N and 6.4%o for <5 18 0. 

To simulate a system alternating between two N z O production 
processes, such as differing relative contributions of nitrification 
and denitrification, with different rates of N 2 0 production and <5 
values, the model was run with both production rate and its 3 
values variable with time (scenarios #4, #5, and #6). The 
production rate and 5 values were adjusted so that the maximum 
rate coincided with the lowest source S values in scenarios #4 and 
7^6 and so that maximum rate coincided with the highest source S 
values in scenario ^5. 

Model Scenario #4: Variable Production Rate, Variable 
Isotopic Composition of Source 

For scenario #4, the resulting N 2 0 concentrations were 
identical to those in model scenario #1, with the maximum 
concentration lagging approximately 2.75 h behind the maximum 
production rate (Figure 7, Table 3). The relationship between the 
5 values of the dissolved and emitted N z O was more complex than 
in other scenarios. The lag time between the maximum source <5 
values and those of dissolved and emitted N 2 0 (when the 
production rate was minimum) was 3.75 h; however, the lag time 
between the minimum source 5 values and those of the dissolved 



and emitted N 2 0 (when the production rate was maximum) was 
only 2.25 h. The difference between the emitted and source N z O 
was 3.1%o to 8.0%o for <5 15 N and 1.3%, to 3.4%o for <5 lfi O. The <5 
values of emitted N 2 0 were closer to those of the source during 
periods of high production rates (and thus higher concentrations) 
than periods of low production rates. However, the flux-weighted 
average 3 values of emitted N 2 0 were equal to the average 
production-weighted source S values because the system was at 
dynamic steady state. 

Model Scenario #5: Variable Production Rate, Variable 
Isotopic Composition of Source 

The isotopic counterpoint to scenario #4 is adjusting the timing 
of maximum N 2 0 production to coincide with the highest S values 
of production (scenario #5). All other parameters were the same 
as scenario #4 (Table 3). The resulting pattern for the S values of 
dissolved N 2 0 was very different than scenario #4 (Figure 8, 
Table 3). While the dissolved N 2 0 concentrations were identical 
to the model scenario #4, the S values of dissolved N 2 0 were 
nearly constant with time. The relationship between the (5 values 
of emitted and source N 2 0 was similar to scenario #4, although 
the instantaneous difference in 3 values were slightly greater. The 
5 values of the dissolved N 2 0 were greatly dampened by the fact 
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Figure 8. Model scenario #5 - Isotopic composition of dissolved and emitted N 2 0 with a variable production rate and variable 
isotopic composition of the source. Maximum production rate is in sync with the greatest <5 15 N and <5 ,8 0 values of the source. 
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that maximum production rate coincided with source 3 values that 
were closest to tropospheric N2O. In scenario #4, the high rates of 
N 2 0 production at 3 values very different than tropospheric N z O 
increased the amplitude of the 3 values of dissolved N 2 0. 

Model Scenario #6: Variable Production Rate, Variable 
Isotopic Composition of Source 

To determine the effects of a lower k on model scenario #4, k 
was reduced from 0.3 m/h from 0.1 m/h for scenario #6. As 
shown above, lower k increased the dissolved N z O concentrations 
and dampened the diel range of 3 values of both dissolved and 
emitted N 2 0 (Table 3, Figure 9). Lower k also increased the lag 
time between the 3 values of emitted and source N 2 0 and 
increased the difference between the 3 values of emitted and 
source N z O (Figure 9). As in all scenarios, the flux-weighted 
average 8 values of emitted N 2 0 were equal to the average 
production-weighted source 3 values. 

Grand River 

The ability of SIDNO to reproduce measured patterns of N 2 0 
concentration and 3 values in a human-impacted river was also 
assessed. The Grand River is a seventh-order, 300 km long river 
that drains 6800 km in southern Ontario, Canada, into Lake 



Erie, see [36,37,51]. There are 30 WWTPs in the catchment and 
their cumulative impact can be observed via the increase in 
artificial sweeteners in the river [52]. 

Samples were collected approximately hourly for 28 h at two 
sites in the central, urbanized portion of the river: sites 9 and 1 1 in 
[51,52]. The upstream site, Bridgeport, is where the river enters 
the urban section of the river at the city of Waterloo and is 
immediately above that city's WWTP. Blair is 26.6 km down- 
stream of Bridgeport and below the cities of Waterloo and 
Kitchener. It is also 5.5 km downstream of the Kitchener WWTP. 
Average river depth at both sites was 30 cm. Values of k were 
determined by best-fit modelling of diel 0 2 and 3 0-0 2 values at 
the sites [36], N 2 0 concentration analyses were performed on a 
Varian CP-3800 gas chromatograph with an electron capture 
detector and isotopic ratio analyses were performed on a GV 
TraceGas pre-concentrator coupled to a GV Isoprime isotope 
ratio mass spectrometer, see [5] for analytical details. 

Data from upstream and downstream of large urban waste- 
water treatment plants on the Grand River show diel patterns in 
N 2 0 saturation and 8 values (Figures 10 and 11). At the 
Bridgeport site, the diel patterns of N 2 0 saturation and (5 15 N 
values were opposite of each other, that is, when N 2 0 saturation 
was highest around sunrise the 3 'N values were lowest and when 
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Figure 9. Model scenario #6 - Isotopic composition of dissolved N 2 0 and emitted N 2 0 with a variable production rate and variable 
isotopic composition of the source. Maximum production rate is in sync with the lowest <5 15 N and <5 18 0 values of the source, k is reduced from 
0.3 m/h to 0.1 m/h. 
doi:1 0.1 371 /journal.pone.0090641 .g009 



when N2O saturation was lowest around before sunset the <5 5 N 
values were greatest. R 2 values between field and model data for 
N 2 0 saturation, <5 15 N, and <5 18 0 values are 0.83, 0.68, and 0.30. 
Model results reproduce the range and sinusoidal patterns of the 
field data though the <5 I8 0 fit was poor in the second half of the 
field data. The diel pattern in S O values was similar to that of 
(5 15 N but was shifted earlier by about 4 h. These patterns were 
similar to those of scenario #4 (variable N2O production and 
variable 3 values of the source N 2 0 coinciding when maximum 
production rates coincided with lowest source 5 values) and the 
result of consistent diel five-fold variability in N 2 0 production and 
variability in the <5 18 0 of the N 2 0 produced in the river. 

At the downstream Blair site, both c5 L, N and <5 18 0 values were 
much lower and exhibited a greater range than at Bridgeport. R 2 
values between field and model data for N 2 0 saturation, 8 N, 
and <5 O values are 0.78, 0.53, and 0.03. Model results reproduce 
the range and peak-and-trough pattern of the N 2 0 saturation and 
f) L, N data. Model results reproduce the range of <5 18 0 values but 
the pattern is not well reproduced. While all data at Bridgeport 
exhibited smooth, sinusoidal diel changes, the data at Blair show 
rapid changes. The diel patterns of N 2 0 saturation and 6 N 
values were opposite of each other, that is, when N 2 0 saturation 



was highest around midnight, the <5 15 N values were lowest and 
when when N 2 0 saturation was lowest during mid-day, the <5 5 N 
values were greatest. The diel pattern in <5 18 0 values was more 
complex at Blair than at Bridgeport suggesting that daytime and 
nighttime were associated with different S O values of N 2 0 
production. These patterns were similar to those of scenario #5 
(variable N 2 0 production and variable <5 values of the source N 2 0 
coinciding when maximum production rates coincided with 
highest source S values) and the result of a five-fold variability in 
day-to-night N 2 0 production and variability in S N and <5 O of 
the N 2 0 produced in the river. 

For both Bridgeport and Blair data, the cause of poorer fits for 
<5 O than S J N deserve further research. Adding concomitant 
measurements of S ''N and 5 O values of NOj~ may provide 
clues about N cycling and help explain some of the observed 
variability in N 2 0 [53]. Predicting <5 18 0-N 2 0 values from 
nitrification [1 1] and denitrification [54] is difficult because of 
the complex relationship between S 0-H 2 0 values and & O- 
N 2 0 values. Additionally, diel variability in N 2 0 reduction to N 2 
[45,46], may also manifest itself in <5 18 0-N 2 0 values because of 
the strong O isotope fractionation factor during denitrification 
[55]. 
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Figure 10. Diel variability in N 2 0 concentration and 5 values at 
Bridgeport in the Grand River, Canada. The time axis begins at 
00:00 on 2007-06-26. Maximum production rate is in sync with the 
greatest <5 18 0 values of the source, while <5 15 N of the source was 
constant. R 2 values between field and model data for N 2 0 saturation, 
<5 15 N, and <5 18 0 values are 0.83, 0.68, and 0.30. This is similar to model 
scenario #4 (Figure 7). 
doi:10.1371/journal.pone.0090641.g010 

Discussion 

Calculating the 8 values of emitted or source N 2 0 is critical for 
regional and global N 2 0 isotopic budgets and also provides 
information about the source of N 2 0 and thus N cycling processes. 
However, SIDNO can simulate the relationships between the 8 
values of dissolved, source, and N z O emitted from aquatic 
ecosystems to the atmosphere. In systems with N 2 0 production 
at dynamic steady state, the S values of dissolved N 2 0 will not 
always be directly indicative of the 8 values of the source N 2 0. 
The difference between dissolved and source 8 values increases as 
N 2 0 saturation decreases (as demonstrated in Figures 2 and 3). 
Even above 1000% saturation (from high production rates and/or 
low k), the 8 values of dissolved N 2 0 will only approach 8 values 
of the source but offset by 0.7%o for <S 15 N and 1.9%o for <5 ls O, a 
result of the a ev values (Figure 3). At constant N 2 0 production 
rates and 8 values, the source and emitted 8 values can be 
quantified since the <5 values of emitted N 2 0 must be identical to 
those of the source and can be calculated from dissolved values 
(Figure 3; equations 5 and 6). 

Our modelling results identified the limitations associated with 
simple interpretation of dissolved N 2 0 isotope data since the 8 
values of dissolved and emitted data are synchronous but rarely 
offset by a constant value. If N 2 0 saturation changes with time, 
the N 2 0 production rate must also have changed with time, 




Time (hours) 

Figure 11. Diel variability in N 2 0 concentration and S values at 
Blair in the Grand River, Canada. The time axis begins at 00:00 on 
2007-06-26. Maximum production rate is in sync with the lowest <5 15 N 
and <5 1S 0 values of the source. R 2 values between field and model data 
for N 2 0 saturation, <3 15 N, and <5 18 0 values are 0.78, 0.53, and 0.03. This is 
similar to model scenario #5 (Figure 8). 
doi:10.1 371/journal.pone.0090641 .g01 1 

provided k had been constant (compare model scenarios # 1 and 
#2 in Figures 4 and 5). In contrast, changes in 8 values of 
dissolved N 2 0 do not require a change in the source 8 values 
(model scenario #1 and Figure 4), while constant 8 values of 
dissolved N 2 0 do not require constant source 8 values (for 
example model scenario #5 in Figure 8). 

When 8 values of the source N 2 0 are variable, the relationship 
between emitted and source N 2 0 becomes complicated. The 8 
values of emitted N 2 0 will lag behind those of the source and the 
amplitude of the diel range of 8 values will be dampened relative 
to the source. The amount of lag and dampening is a function of k, 
N 2 0 production rate and timing, and the proximity of the source 
8 values to those of the atmosphere (compare Figures 2 with 3 and 
Figures 4 with 6). Qualitatively, the 8 values of emitted N 2 0 will 
be similar to the source if the equilibration time of dissolved N 2 0 
is small relative to the period of source variability (e.g., 24 h period 
due to diel changes in N cycling [36,45]). Assuming homogeneous 
N 2 0 release upstream, the equilibration time can be approximat- 
3z 

ed from a decay curve as — , where z is mean depth [49] . If z is 
K 

small and/ or k is high, the equilibration time will be short and the 
8 values of the emitted N 2 0 will be close to the source. With 
decreasing k (or increasing equilibration time), the 8 values of 
emitted N 2 0 will lag farther behind and will always have a smaller 
range of 8 values than the source. At the most extreme case, the 
variability in the 8 values of emitted N 2 0 will be reduced to nearly 
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Table 4. Summary of the results of the SIDNO modelling as the predictive relationship between observations and implications. 



Observed Parameter 


Implications 


Examples 


Dissolved N 2 0 concentration is constant with time 


The N 2 0 production rate is constant with time (if k 
and temperature are also constant). The N 2 0 flux to 
the atmosphere is equal to the production rate. This 
may not be true if the concentration is close to 
atmospheric equilibrium. 


Scenarios #2 and #3 


Dissolved N 2 0 concentration is variable with time in 
a sinusoidal pattern 


The N 2 0 production rate is variable with time (if 
concentration change cannot be explained by 
change in k or temperature). The average N 2 0 
flux to the atmosphere is equal to the average 
production rate. 


Scenario #1 


<5 15 N and <5 18 0 of dissolved N 2 0 is constant with time 


The observation is inconclusive. At concentrations 
near atmospheric equilibrium, isotopic composition 
of dissolved N 2 0 will approximate tropospheric N 2 0, 
egardless of source values. A constant isotopic signature 
of dissolved N 2 0 that is different from tropospheric N 2 0 
can indicate either a constant source (if production rate 
is constant), or a variable source. 


Scenario #5 


<5 15 N and <5 18 0 of dissolved N 2 0 is variable with time 
(slope of data on a <5 18 0— c) 15 N cross-plot tends toward 
tropospheric N 2 0} 


The change in <5 15 N and <5 18 0 of dissolved N 2 0 is likely a 
result of a change in concentration, but it is possible that 
the source is variable with time, if the t) 15 N and <5 18 0 values 
of the source also trend through the value for tropospheric N 2 0. 


Scenarios #1, #2, 
^f4, ana fpo 


Calculated f3 15 N and t) 18 0 of emitted N 2 0 is constant 
with time 


The isotopic composition of the source is constant with time, 
and equal to the calculated value for emitted N 2 0 


Scenario ^1 


Calculated c5 15 N and <5 18 0 of emitted N 2 0 is variable 
with time 


The isotopic composition of the source is variable with time. 
The range in <5 15 N and 6 ]8 0 of emitted N 2 0 is the minimum 
for the range in that of the source. The flux weighted 
average c5 15 N and t> 18 0 of emitted N 2 0 is equal to the 
production weighted average source values. 


Scenarios ^2—^6 


Long residence time relative to variability of source 
(need to independently determine k) 


The changes in N 2 0 concentration and f3 15 N and t> 18 0 of 
emitted N 2 0 will be dampened relative to, and lag behind, 
that of the source 


Scenarios #3 and #6 


Short residence time relative to variability of source 
(need to independently determine k) 


The changes in N 2 0 concentration and (5 15 N and t> 18 0 of 
emitted N 2 0 will be indicative of changes in the source 


Scenarios #1, #2, #3 
and #4 



doi:1 0.1 371 /journal.pone.0090641 .t004 



zero and <5 values of the emitted N 2 0 would be equal to the 
average production-weighted source S values. At very long 
equilibration times, the probability of N 2 0 consumption increases, 
a process not explicitly included in SIDNO where the S value of 
the source N 2 0 is simply that which is released to the water 
column. 

Separating N z O production into nitrification and denitrification 
requires independent knowledge about the S values of the source 
N and O in aquatic ecosystems. It is therefore not possible to state 
a single c5 15 N value for nitrification-N 2 0 and one for denitrifica- 
tion-N 2 0 applicable to all aquatic ecosystems. The <5 1S N value of 
the N 2 0 precursors NH^ and NO^ vary across ecosystem as a 
result of human impact and N loading (agricultural and WWTP) 
as well as the source of N, and additional N transformations in the 
aquatic ecosystem. For example, along the length of the Grand 
River, 5 5 N values of NH^ and NOj~ exhibit systematic trends 
resulting from the confluence of agricultural tributaries and large 
urban waste-water treatment plants (Scruff et al, unpublished 
results, [53]). Nevertheless, these values can be measured and 
biogeochemical relationships between N species, redox, and N 2 0 
can be used as supporting information for process separation (e.g., 
[5,12,36,45]). The 5 18 0 value of N 2 0 will also vary across 
ecosystems as a result of its close relationship with <5 18 0-H 2 0 and 
to a lesser extent c5 1!i O-0 2 [10,11,54,56]. Fortunately, (5 18 0-H 2 0 
values can be easily predicted and measured [57]. Thus, once <5 
values of N 2 0 precursors have been identified, biogechemical data 
can provide an indication about the diel pattern of N 2 0 



production processes, and ranges of potential end-member 5 
values can be calculated (e.g., [5] summarize isotopic fractionation 
for <5 15 N and [10,1 1,54,56] for <5 1!i O) and the model used to fit the 
field data. 

Conclusions 

In aquatic ecosystems, the instantaneous 3 values of N 2 0 
emitted to the atmosphere are easily calculated if the water 
temperature and dissolved N 2 0 concentration and 8 values are 
known. Our modelling efforts illustrate that complex relationships 
exist between dissolved and source N 2 0 and that the S values of 
dissolved N 2 0 are not always representative of either the N 2 0 
produced or emitted to the atmosphere. Thus, calculated d values 
of the emitted N 2 0 are the values that should be used to draw 
conclusions about N 2 0 emission from aquatic systems and the 
global N 2 0 cycle rather than the more commonly used 
instantaneous values (Table 4). The flux-weighted (5 values of 
emitted N 2 0 can provide average production-weighted c5 values 
of the N 2 0 source under dynamic steady-state in aquatic 
ecosystems. 

If the 8 values of emitted N 2 0 are constant with time, either the 
<5 values of the source must also be constant or the N 2 0 
equilibration time is very long. However, if the calculated S values 
of emitted N 2 0 vary with time then the <5 values of the source must 
also vary with time producing a diagnostic pattern. These findings 
are more robust than using dissolved S values alone since dissolved 
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S values can change simply with a change in N 2 0 production rate, 
changes in source S values, and changes in k. N 2 0 residence time, 
dependent on production rate, k, and z, will determine the lag 
time between the S values of emitted and source N 2 0. The 
difference in timing between maxima and minima S values of 
emitted N 2 0 and the maxima and minima of dissolved N 2 0 is 
indicative of how the S values of the source change. For all these 
reasons, we urge caution when using single samples of N 2 0 
concentration and S to calculate fluxes of N 2 0 to the atmosphere 
and inferring N 2 0 production pathways. 

Ultimately, the dynamic model SIDNO may be used to estimate 
N 2 0 production rate and S values of the N 2 0 source, an 
indication of the production pathway and N cycling, in aquatic 
ecosystems via inverse modelling. If physical properties, such as 
depth and temperature are known, SIDNO may be used to fit the 
measured field data (N 2 0 concentration and S values) by adjusting 
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